Gravitational waves from black hole collisions via an eclectic approach 
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We present the first results in a new program intended to make the best use of all available technolo- 
gies to provide an effective understanding of waves from inspiralling black hole binaries in time for 
imminent observations. In particular, we address the problem of combining the close-limit approxi- 
mation describing ringing black holes and full numerical relativity, required for essentially nonlinear 
interactions. We demonstrate the effectiveness of our approach using general methods for a model 
problem, the head-on collision of black holes. Our method allows a more direct physical understand- 
ing of these collisions indicating clearly when non-linear methods are important. The success of this 
method supports our expectation that this unified approach will be able to provide astrophysically 
relevant results for black hole binaries in time to assist gravitational wave observations. 
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Binary black hole systems pose one of the most excit- 
ing and challenging problems of general relativity, con- 
stituting not only a fundamental gravitational two-body 
problem, but also an important astrophysical problem 
of direct and immediate observational relevance. Grav- 
itational waves from binary black hole mergers are con- 
sidered one of the most promising candidates for experi- 
mental detection by the first wave of large interferomet- 
ric gravitational wave observatories coming on line over 
the next few years. These imminent observations present 
an urgent call to the theoretical relativity community to 
immediately provide any information possible about the 
radiation that might be expected from these collisions. 

The problem divides physically into three phases. Ini- 
tially, a slow adiabatic inspiral lasting nearly until the 
black holes are so close the orbital motion destabilizes, 
a brief period of strong essentially non-linear two-body 
interaction, and the linear ring-down of the newly formed 
remnant black hole to stationarity. Correspondingly, the- 
orists have approached the problem along three primary 
avenues: the post-Newtonian (PN) slow-motion approx- 
imation applicable in the inspiral phase, the 'close-limit' 
(CL) single perturbed black hole approximation handling 
the ring-down, and the full numerical simulation (FN) of 
Einstein's equations, which could ideally handle the en- 
tire problem on a large computer, but is so-far limited to 
brief evolutions on small 3D domains. Nevertheless the 
full numerical approach should be vital to treating the 
intermediate, essentially non-linear phase. 

In order to form the best theoretical model possible 
for the radiation from these systems we feel it is vital to 
combine these three approaches focusing the numerical 
simulations squarely on the intermediate phase of the in- 
teraction where no perturbative approach is applicable. 
The state-of-the-art in these three fields has advanced to 
the point where we can expect such an eclectic approach 
to provide a reasonable model for binary black hole ra- 
diation without depending on further advancements. In 



the cases where it has been applied the close-limit model 
has proven to be a reliable model for radiation after the 
system has formed a common event horizon, and work 
in this field has advanced sufficiently so that arbitrary 
perturbations can now be calculated routinely M . In full 
numerical relativity, parts of the plunge of rather gen- 
eral black hole systems, the grazing collision of two black 
holes with linear momenta and spins, can be simulated 
H . And the post-Newtonian method has advanced to the 
point where it might be trusted even for black holes ap- 
proaching the last stable orbit (LSO). Recent estimates 
H suggest that in the absence of spins there are 0.6 or- 
bits left for full- numerical treatment, and this part of the 
plunge (roughly 50M evolution time), numerical relativ- 
ity should be able to handle today. The primary obstacle 
to proceeding with the combined model is the construc- 
tion of appropriate interfaces between the three exist- 
ing models. A recent interest within the post-Newtonian 
and gravitational wave research community in providing 
Cauchy data for simulations may soon solve the problem 
of the PN-FN interface J|,^| . In this letter we introduce 
a general approach to providing the FN-CL interface. 

The nominal result of a numerical simulation of Ein- 
stein's equations is a time succession of values on a 3D 
grid for the spatial metric and extrinsic curvature hold- 
ing the geometric spacetime information. For binary 
black hole simulations we expect the late time behav- 
ior of the system to be best characterized as a "ring- 
ing" black hole with outgoing radiation, with perturba- 
tion theory providing a good model for the dynamics. 
The perturbative model not only allows an inexpensive 
continuation of the evolution, but also supplies a clear in- 
terpretation of the dynamics not manifest in the generic 
numerical simulation. The dynamics reduces to the evo- 
lution of a single complex field, the Newman-Penrose 
Weyl scalar ip± = C a (i^&n a fhP rf m s , obeying a linear hy- 
perbolic equation. Because of the axisymmetry of the 
background Kerr black hole the problem can be further 
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simplified by Fourier decomposition of tpi, reducing to a 
series of 2D evolution problems for the axial mode compo- 
nents of ip4 evolving according to the Teukolsky equation 
1- 

Several important steps are required to concretely 
implement the FN-CL interface: 1) Specify the back- 
ground black hole by its mass M and angular momen- 
tum a = J I'M 2 . 2) Construct a space-like slice from the 
late-time region of the numerical spacetime which will 
be mapped to a constant time slice in the perturbative 
calculation. In general this slice may not be related to 
the numerical foliation. 3) Specify an embedding explic- 
itly mapping the numerical slice to the corresponding 
slice in the background spacetime. 4) Specify a (null and 
complex) tetrad, (i M , n M , m M , fh^) , on the numerical slice 
which will map, on the background slice, to an approx- 
imation of the standard tetrad used in the perturbative 
formalism where and are conveniently chosen to 
lie along the two-degenerated principal null directions of 
the background spacetime. 5) Using the specified tetrad 
and the numerical data for the metric gij and the ex- 
trinsic curvature 2£y on the slice calculate ip^ and dtip4- 
These will provide the Cauchy data for the perturbative 
evolution. 6) Evolve with the time-domain Teukolsky 
equation to determine the subsequent perturbative dy- 
namics. The heart of the problem is making the specifi- 
cations required in (2-4). The general idea is to numer- 
ically compute physical quantities or geometrical invari- 
ants and relate them to their analytic expressions in the 
perturbatively preferred coordinate system. Coordinate 
information can also be computed, for instance, dragging 
information along geometrically defined trajectories (like 
geodesies) from less problematic 'faraway' parts of the 
spacetime where the two slices are much closer to each 
other. Another possibility is to use the local null struc- 
ture of the spacetime to determine the eigenvalues of the 
Weyl tensor, perturbatively related to the principal null 
directions of the Kerr background. There is generally 
no unique way to make these specifications but the first 
order gauge and tetrad invariancc of the perturbative for- 
malism implies that the results will not depend strongly 
on small variations in these choices. Step (5) was explic- 
itly worked out in Ref. H . 

Since a concrete implementation requires us to make 
choices for which there is no clear mathematical prefer- 
ence, we will proceed by trying first the simplest pos- 
sible specifications and adding sophistication only when 
it seems to be necessary. We begin with a model bi- 
nary black hole problem which has already been solved 
by 2D numerical relativity and close-limit perturbation 
theory, head-on collisions of initially resting equal-mass 
black holes (Misner initial data). At the same time, we 
will try not to tune our techniques too closely to this par- 
ticular example so that our method can be readily gen- 
eralized. For this reason we will perform our numerical 
evolutions generically in 3D, using well-tested, numerical 



techniques and codes (Cactus M) "off the shelf" with no 
fine-tuning for this problem. We also apply perturbation 
theory as described by the Teukolsky equation, allowing 
for a rotating black hole background, without multipole 
decomposition. Specifically, for the numerical evolutions 
we have used the ADM system of Einstein's equations 
with maximal slicing for the lapse and vanishing shift, fi- 
nite differenced on a 128 3 (octant mode) numerical grid, 
initially mapped non-uniformly to the standard Misner 
coordinates to allow a distant outer boundary without 
sacrificing resolution in the inner region. We express the 
Teukolsky equation in Boyer-Lindquist coordinates, al- 
though it may be convenient in the future to evolve the 
perturbations in another gauge, such as a Kerr-Schild 
representation of the Kerr metric || . 

We implement the steps listed above as follows: 1) In 
this case there is no angular momentum so the back- 
ground reduces to Schwarzschild, a = 0. Since only 
about 0.1% of the system's mass will be lost as radi- 
ation we specify the background mass as equal to the 
initial ADM mass. 2) We make the simplest choice of 
background slice by identifying the numerical slice with 
a Schwarzschild time slice. Numerical experience with 
Schwarzschild black hole evolutions in maximal slicing 
suggest a strong correspondence. 3) For the embedding, 
it is clear that the trivial choice, identification of numer- 
ical and background coordinates is inadequate because 
the black hole horizon must invariably expand in this 
numerical gauge. On the other hand the same expan- 
sion has the tendency to drive the exterior region toward 
manifest spherical symmetry. A reasonable estimate for 
the map into the background Schwarzschild slice is a triv- 
ial identification of the numerical and background 9 and 
(j) coordinates, with some relabeling of the constant-r 
spheres. We account for the radial rescaling by choos- 
ing the background radius r' so that the value of Weyl- 
curvature invariant X = C a bcdC ahcd averaged over in 
the numerical slice coincides with its background value 
I = 3M 2 /r' e in the background slice. 4) We define an 
appropriate, manifestly orthonormal, tetrad primarily by 
identifying timelike normal, radial, and azimuthal direc- 
tions. The unit normal and radial direction vectors pro- 
viding the spatial components of P and n M . The complex 
vectors and fh} 1 point within the spherical 2-surface. 
At each step, a Gram-Schmidt procedure is first used 
to ensure that the triad remains orthonormal. Then a 
type III (boost) null rotation fixes the relative normaliza- 
tion of the two real- valued vectors to make it consistent 
with the tetrad assumed in the perturbative calculation. 
5) Within the full numerical simulation we compute the 
values ipi an d dtip4 consistent with our tetrad specifica- 
tion using the formulas in Ref. |(| and interpolate these 
Cauchy data (using splines) to generate data directly us- 
able by the Teukolsky code developed in Ref. || . For the 
perturbative evolutions we use —20 < r*/M < 50 with 
n g x n r * = 48 x 700. 
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FIG. 1. The total radiated energy from two black holes 
initially at rest (Misner initial data). The solid line repre- 
sents Price-Pullin results Q (labeled as PP, £ = 2) for the 
close limit approximation. 2D full numerical results [[[o| arc 
here given by full circles and the new FN-CL results are la- 
beled with stars. The minimal amount of evolution needed for 
fio = 1.8, 2.2, 2.7 (L « 3.4, 4.5, 6.4) to reach the perturbative 
regime are approximately 1M, 3M, 6M respectively. 



We evolved Einstein's equations numerically from Mis- 
ner initial data for several different initial separations la- 
beled by the parameter /io (corresponding to proper sep- 
arations L as shown in Fig. |l|) . A typical duration of the 
total evolution was t = 10M and we extracted Cauchy 
data every t = 1M. A transition time tx was deter- 
mined by methods detailed below. After each Teukolsky 
code evolution we extract the full relevant signal of the 
waveforms, which typically lasted for t < lOOAf. The 
resulting radiation energies are shown in Fig. [j] where 
we compare our 3D results with the results of Ref. ]Tc| ] 
where explicit use of the symmetries of the problem have 
been implemented in a 2D simulation. The other case for 
comparison is the Price-Pullin curve providing the 
pure close- limit result. While all three predictions agree 
very well for small initial proper separations L/M < 3, 
it is clear that for larger separations the close limit and 
full numerical curves deviate considerably. Our results 
follow quite precisely the 2D computations. A minimal 
full numerical evolution time (given by our linearization 
time below) is essential in obtaining the above agreement. 
Evolution of exact initial data only perturbatively does 
not reproduce the full numerical results for large separa- 
tions, but follows the PP curve [G3. 

Extracting waveforms every IM of non-linear numer- 
ical evolution allows us to study the transition to linear 
dynamics, and to perform important consistency tests 
on our results. If we have made a good definition of the 
perturbative background in steps (1-4) above then we 
can expect our radiation waveform results to be indepen- 
dent of the transition time, tx, once the linear regime is 
reached and for as long as the numerical simulation COn- 
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FIG. 2. Radiated energy versus transition time. These fig- 
ure show a clear plateau after linearization until numerical 
error begins to cause problems after 6M. 



tinues to be accurate. We apply two independent criteria 
for estimating the onset of linear dynamics, the speciality 
invariant prediction based only on the Cauchy data and 
another estimate based on the stability of the radiation 
waveform phase. The speciality invariant introduced in 
]l3[ predicts linear dynamics when S — 21J 2 jT? differs 
from its background value of unity by less than a factor 
of two outside the (background) horizon. Such a devia- 
tion from algebraic speciality implies significant "second 
order" perturbations. The phase of the radiation also 
provides an indicator of linear dynamics. Starting with 
detached black holes, we expect an initial period of weak 
bremsstrahlung radiation followed by the appearance of 
quasinormal ringing. On the other hand, switching to 
perturbative evolution prematurely leads to immediate 
ringing. Hence we first observe a series of phase delays 
for the beginning of the ringing until the actual ringing 
takes place, thereafter no phase shift should be seen. The 
value of tx when the phase freezes gives a precise esti- 
mate of time for linearization of the system. We find that 
both estimates for linearization time are in good agree- 
ment, yielding that t < 1M for the /j,q = 1.2 case, t « 1M 
for = 1-8, 3M for /i = 2.2, and 6M for [i = 2.7. The 
linearization time is somewhat longer than the "ringing 
times" reported in Fig. 7 of Ref. jl0| indicating that lin- 
earization occurs slightly after the onset of "ringing" for 
the stronger collisions. Our linearization times are still 
much shorter than those for the appearance of a common 
apparent horizon. 

Two example curves of energy versus tx are shown in 
Fig. g. Before the linearization time the premature appli- 
cation of perturbation theory tends to result in an overes- 
timate of the energy. After linearization there is a plateau 
region when the energy is insensitive to tx as is required 
if we have defined a useful FN-CL interface. Eventually, 
after 6M in these cases, numerical errors caused by the 
"grid-stretching" inherent in the use of maximal slicing 
tend to result again in an overestimate of the energy. A 
stronger indication of the robustness of our method is ev- 
ident in the waveforms themselves. After linearization, 
the waveforms should also be independent of tx- Fig. || 
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FIG. 3. Waveforms and Cauchy data for initial separation 
fio — 2.2. The observer location is r* = 20M and 9 = 7r/2. 
The data have been extracted after t = 4M and 5M of full 
3D nonlinear evolution respectively. There are evident differ- 
ences in the extracted Cauchy data at the two times because 
of the intervening evolution. Nevertheless the resulting wave- 
forms agree, demonstrating the equivalence of the linear and 
non-linear evolution through this period, as well as the ro- 
bustness of our methods. 



shows an example of this comparing tr = 4M and 5M 
for the no = 2.2 case. Despite the fact that the Cauchy 
data at transition time is very different, the waveform is 
almost identical. The waveform quite agrees (apart from 
the reversed sign) with the ^4 published in Ref. j[o|, Fig. 
13. It is worth noting here that our waveforms for the 
Misner data seem to be the first complete ones computed 
using 3D full numerical relativity. 

Perturbation theory is very useful to gain information 
about waveforms from numerical spacetimes. Customar- 
ily this interface has been implemented only on a time- 
like surface to determine radiation content "far-away" 
from the black hole. A much more natural boundary 
between the linear and non-linear regimes occurs on a 
spacelike interface defined by the time beyond which non- 
linear black hole perturbations no longer contribute sig- 
nificantly to the radiation. We have taken a general ap- 
proach to the problem of providing such a FN-CL in- 
terface which we believe is essential to providing timely 
estimates of binary black hole waveforms. We are aware 
of only one previous attempt to make a combined use of 



numerical and close-limit evolution implemented in the 
case of two black holes formed by collapsing matter [jlJJ , 
using a 2D code and I — 2 metric perturbations (a la Zcr- 
illi) of the Schwarzschild background. Our method aims 
toward complete generality using full 3D numerical sim- 
ulations and applying perturbation theory as described 
by the Teukolsky equation, applicable to arbitrary rem- 
nant black hole backgrounds. This approach is directly 
applicable to a unified eclectic model of colliding black 
holes joining the close-limit, full numerical relativity, and 
post-Newtonian methods. To our knowledge this is the 
first time such an approach is proposed and turned into 
a concrete and generic scheme. The success in this test 
case encourages our hopes for a providing theoretical re- 
sults on black hole merger waveforms in time to assist the 
first gravitational wave observations. We will direct our 
future work toward a fully combined PN-FN-CL model 
for estimating astrophysically relevant binary black hole 
collision waveforms. 
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